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Abstract 

We consider the stability of synchronized states (including equilibrium point, 
periodic orbit or chaotic attractor) in arbitrarily coupled dynamical systems 
(maps or ordinary differential equations). We develop a general approach, 
based on the master stability function and Gershgorin disc theory, to yield 
constraints on the coupling strengths to ensure the stability of synchronized 
dynamics. Systems with specific coupling schemes are used as examples to 
illustrate our general method. 
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Large networks of coupled dynamical systems that exhibit synchronized static, periodic 
or chaotic dynamics are subjects of great interest in a variety of fields ranging from biology 
[1] to semiconductor lasers [2] to electronic circuits [3]. For a given problem it is essential 
to know the extent to which the coupling strengths can be varied so that the synchronized 
state remains stable. Early attempts [4,5] at this question have typically looked either at 
systems of very small size or at very specific coupling schemes (diffusive coupling, global all 
to all coupling etc. with a single coupling strength). Recent work [6,7] introduced the notion 
of a master stability function that enables the analysis of general coupling topologies. This 
function defines a region of stability in terms of the eigenvalues of the coupling matrix. In 
this paper we present a general method that provides explicit constraints on the coupling 
strengths themselves by combining the master stability function with the Gershgorin disc 
theory. Our approach is applicable to both coupled maps and coupled ordinary differential 
equations (ODEs). Commonly studied coupling schemes are used as illustrative examples. 

Coupled Maps. The system we consider is represented by 

1 N 

x> + 1) = f (x») + - £ G y • H(tf(n)), (1) 

j=i 

where x l (n) is the M-dimensional state vector of the ith map at time n and H : R M — > R M 
is the coupling function. We define G = [Gy as the coupling matrix where Gij gives the 
coupling strength from map j to map i. The condition y^Gj, = is imposed to ensure that 

j 

synchronized dynamics is a solution to Eq. (1). 

Linearizing Eq.(l) around the synchronized state x(n), which evolves according to x(n + 
1) = f(x(n)), we have 

1 N 

z\n + 1) = J(x(n)) • z»(n) + - £ G i3 • £>H(x(n)) • z», (2) 

JV 3=1 

where z*(n) denotes the ith map's deviations from x(n), J(-) is the M x M Jacobian matrix 
for f and -DH(-) is the Jacobian of the coupling function H. In terms of the M x N matrix 
S(n) = (z 1 (n) z 2 (n) • • • z N (n)), Eq. (2) can be recast as 

S(n + 1) = J(x(n)) • S(n) + ^H(x(n)) • S(n) • G T , (3) 

According to the theory of Jordan canonical forms, the stability of Eq. (3) is determined 
by the eigenvalue A of G. Denote the corresponding eigenvector by e and let u(n) = S(n)e. 
Then 

u(n + 1) = (J(x(n)) + -^A • DH(x(n))) • u(n). (4) 

So the stability problem originally formulated in the M x N space has been reduced to a 
problem in a M x M space where it is often the case that M <g N. It is worth mentioning 
that this eigenvalue based analysis is valid even if the coupling matrix G is defective [8]. 
We note that A = is always an eigenvalue of G and its corresponding eigenvector is 

N 

(11 ■ ■ • 1) T due to the synchronization constraint X) Gij — 0- 111 this case, Eq. (4) can 

3=1 
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be used to generate the Lyapunov exponents for the individual system, which we denote 
by hi = h max > h 2 > • • • > hu- These exponents describe the dynamics within the 
synchronization manifold defined by x* = x Vi. 

The subspace spanned by the remaining eigenvectors is transversal to the synchronization 
manifold, the dynamics in which will be stable if the transversal Lyapunov exponents are 
all negative. To examine this problem, we treat A in Eq. (4) as a complex parameter 
and calculate the maximum Lyapunov exponent [i m ax as a function of A. This function 
is referred to as the master stability function by Pecora and Carroll [6]. The region in the 
(Re(A), Im(A)) plane where fx max < defines a stability zone denoted by f2. Figure 1 shows a 
schematic of two possible configurations of f2. Whether f2 is an unbounded area [Fig. 1(a)] or 
a bounded one [Fig. 1(b)] is contingent on the coupling scheme and other system parameters. 
The origin, which is the zero eigenvalue of G, may or may not lie in the stability zone. For 
example, for equilibrium or periodic state in coupled maps, the origin is in Q, but for chaos, 
it lies outside of Q. We note that, typically, O is obtained numerically. In some instances 
analytical results are possible (see below). 

Clearly, if all the transversal eigenvalues of G lie within Q, then the synchronized state is 
stable. Here we seek constraints applicable directly to the coupling strengths. This problem 
is dealt with by combining the master stability function with the Gershgorin disc theory. 

The Gershgorin disc theorem [9] states that all the eigenvalues of a n x n matrix A = [a^] 
are located in the union of n discs (called Gershgorin discs) where each disc is given by 

{z e C : \z - a u \ < Ki|}> i = l,2,...,n. (5) 

To apply this theorem to the transversal eigenvalues we need to remove A = 0. We appeal 
to a order reduction technique in matrix theory [10] which leads to a (N — 1) x (N — 1) 
matrix D whose eigenvalues are the same as the eigenvalues of G except for A = 0. 

Suppose that, for a given matrix G, we have knowledge of one of its eigenvalues A and 
the eigenvector e. Through proper normalization we can make any component of e equal 
one. Here, without loss of generality, we assume that the first component is made equal 1, 
namely, e = (1, e N _ 1 ) T . Rewrite G in the following block form: 



G = " n ' (6) 




with r = (G 12 , • • • , G 1N ) T , s = (G 21 , • • • , G N1 ) T and 

( G22 ■ • ■ G 2 n ^ 
Gjv_i = : : : 

\ Gn2 ■ ■ ■ Gnn J 



(7) 



Choose a matrix P in the form 

P=( ' r" )■ 




Here Iat-i is the (N — 1) x (N — 1) identity matrix. Similarity transformation of G by P 
yields 
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Since P _1 GP and G have identical eigenvalue spectra, the (N — 1) x (JV — 1) matrix 

D 1 = Ga?_i — eAr_ir T (10) 

assumes the eigenvalues of G sans A. We can obtain N different versions of the reduced 
matrix, which we denote by D fc (k = 1, 2, • • • , N), depending on which component of e is 
made equal 1. 

Applying the above technique to the coupling matrix G by letting A = and e = 
(11 ■ • • 1) T we get D fc = [dij] where d^ = Gij — Gkj- From the Gershgorin theorem the 
stability conditions of the synchronized dynamics are expressed as 

1. The center of every Gershgorin disc of D fc lies inside the stability zone Q. That is, 
(Ga — Gki, 0) G f2. 

2. The radius of every Gershgorin disc of D h satisfies the inequality ^ ~ G ki \ < 

j=l,j¥=i 

S(Gu — G^), i = 1, 2, . . . , N and i ^ k. Here S(x) is the distance from point x on 
the real axis to the boundary of the stability zone Q. 

As k varies from 1 to N, we obtain N sets of stability conditions. Each one provides sufficient 
conditions constraining the coupling strengths. 

Coupled ODE's. The above procedure for obtaining stability bounds can also be applied 
to coupled identical ODEs written as 

1 N 

x* = F(x*) + -£^H(x,), 

iV 3=1 

where x l is the M-dimensional vector of the ith node. The dynamics of the individual node 
is x = F(x). Linearizing around the synchronized state we get 

. l N 

z* = J(x).z' + -£Gy-£>H(x).z», (12) 

iV 3=1 

where z* denotes deviations from x, J(-) and -DH(-) are the M x M Jacobian matrices for 
the functions of F and H. Adopting Jordan canonical form, we obtain 



u = 



J(x) + 1a • DH(x) 



(13) 



where A is an eigenvalue of G. Performing the same analysis as for coupled maps, we obtain 
the same stability conditions as given above. 

Examples. We now illustrate the general approach by applying the above results to 
two examples where analytical results are possible. In the first example we consider the 
coupled differential equation systems with H(x) = x [7]. It is easy to see that DH is a 
M x M identity matrix. The Lyapunov exponents for Eq. (13) are easily calculated since 
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the identity matrix commutes with J(x). Denoting them by /ii(A), /12(A), . . . , /im(A), we 
have 

Vi(\) = h i + -^Re(\), i = l,2,...,M. (14) 

For stability, we require the transversal Lyapunov exponents (A 7^ 0) to be negative. This 
is equivalent to the statement that the maximum Lyapunov exponent is less than zero: 

1 

AWr(A) = h max + ^Re(A) < 0. (15) 

In other words, the stability zone f2 is the region defined by Re(A) < —Nh max . The distance 
function from the center of each Gershgorin disc to the stability boundary is given by: 
5(Gu — G ki ) = —h max — (Gu — G ki ) (i — 1, . . . , TV, i ^ k). Thus the fc-th set of stability 
conditions are 

{Gu - G ki ) < -Nh max , (16) 

N 

E \Gji — Gki\ < —Nh max — (Gu — G^), (17) 

% = 1,2, ...,N, i^k. 

It is obvious that the second inequality implies the first one. So the stability condition for 
the synchronized state (whether an equilibrium, periodic or chaotic state) is given by 

N 

E \G ji -G ki \ + (G ii -G ki )<-Nh max , i = l,2,...,N, i^k. (18) 

j=l,j¥=i 

When the coupling is symmetric, i.e. Gij = Gji, Rangarajan and Ding [11], based on the 
use of Hermitian and positive semidefinite matrices, derived a very simple stability constraint 

G^ > h max , (19) 

We show here that Eq. (19) is a consequence of the more general stability conditions given in 

N 

Eq. (18). This can be seen as follows. First consider k — 1. Substituting Gu = — E Gji 
(synchronization condition) and simplifying we get the following equation: 

N N 

E \Gji-G u \- E G ji -2G li <-Nh max: i^l. (20) 

If Gji — G\i is positive for all allowed % and j values, it is easy to see that the above stability 
condition is satisfied given the condition in Eq. (19). However, if more than two such terms 
are negative we have a problem. We can get around this by considering the other (N — 1) 
sets of stability conditions obtained by setting k = 2, 3, . . . , iV in Eq. (18): 
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N N 

\Gji — G 2 i\ — ^ Gji — 2G 2 i < —Nh max , i ^ 2 

; (2i) 

N-1 N-1 

— G]yi\ — ^ Gjj — 2G]yi < —Nh max , i ^ N 

If we take the average of the inequalities over k, cancellation takes place, resulting in a 
simplified inequality that will be satisfied if the sufficient condition given in Eq. (19) is met. 
In other words, the previously derived stability condition is obtained as a special case when 
we require the coupling strengths to meet the N stability conditions simultaneously. 

In the second example, we consider a coupled map with H = f [5] . Under this assumption, 
DH = J and the linearized equation [cf. Eq. (4)] reduces to 

u(n+l) = (A/iV + l)J(x(n))u(n). (22) 

The Lyapunov exponents for Eq. (22) are easily calculated analytically. Denoting them by 
Ati (A), /i 2 (A), • • • , //m(A), we have 

m(\) = hi + hi\\/N + l\, i = l,2,...,M. (23) 

For stability, we require n max (X) = h max + ln \\/N + 1| < 0. In other words, the stability 
zone is defined by 

\\ + N\ < N exp{-h max ) (24) 

The distance from the center of each Gershgorin disc to the boundary is easily calculated 
to be S(Gii - G ki ) = Nexp(-h rnax ) - \N + G u - G ki \ (i = 1, . . . , N, i ± k). Thus the 
conditions of stability are 

N 

E - + \N+ (Gu - G ki )\ < Nexp(-h max ), (25) 
i — 1, . . . , TV, % ^ k, k = 1 or 2 or • • • or N. 

For each k from 1 to N, we obtain a set of sufficient stability conditions. 

In [11], a simple stability bound for synchronized chaos in the case of symmetric coupling 
was obtained as 

[1 - exp(-/i ma:E )] < Gij < [1 + exp(-/wr)], Vi, j. (26) 

This can again be derived from the general stability condition in Eq.(25) with the averaging 
technique used above. 

In summary, we have set up a general formalism to study the stability of synchronized 
state in coupled identical maps and ordinary differential equations. We have also considered 
the often used coupling function for coupled maps and coupled ODEs and given analytical 
results in these cases. We have also shown that known stability bounds can be derived from 
our more general results. 

The work was supported by US ONR Grant N00014-99-1-0062. GR was also supported 
by the Homi Bhabha Fellowship. 
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FIGURE CAPTION 

Figure 1: Schematic illustrations of the stability zone, (a) unbounded area, (b) bounded 
area. 
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